close all;
clear all;
l0=0.25; % rest length m
m = 0.1; % 1 kg mass
k = 1500; %spring stiffness
b = 10;%m/seconds
tMax = 0.5; %seconds
deltaT = 0.001; %seconds
FappTot = 0.5; %Newtons
nMax = [2];

%x1 = l0 ; % initial extention
%xDot1 =0; %initial velocity
%xDoubleDot1 = 0;

%x2= 2*l0 ; % initial extention
%xDot2 =0; %initial velocity
%xDoubleDot2 = 0;

x = zeros(nMax,1);

for(i=1:nMax)
	x(i) = i*l0;
endfor
disp(x);
xDot = zeros(nMax,1)
xDoubleDot = zeros(nMax,1)

store = zeros((tMax/deltaT)+1, nMax);
Fapp = zeros(nMax,1)

t = 0;
%Fapp2 = Fapp;
%Fapp = 0;
Fctr = 1;
Fapp(1) = FappTot
disp("Time loop starting");
store1 = x(1);
store2 = x(2);

while (t<tMax)

	%disp(t);		
%	if (( t >(tMax) && (Fapp[Fctr] != Fapp))
%		Fapp1=Fapp;
%		disp("Fapp1 has changed");
%	endif
%	disp("starting intergration");
	%	xDoubleDot1 = (1/m)*( - Fapp1 - ( k*(x1-l0)) - (b*xDot1) + k*(x2-(x1+l0)) );
	xDoubleDot(1) = (1/m) * ( -Fapp(1) - (k*(x(1)-l0)) - (b*xDot(1))  + k*(x(2) - (x(1) +l0)) );
	%Integration of Positionn and velocity;
	xDot(1) = xDot(1) + xDoubleDot(1)*deltaT;
	x(1) = x(1) + xDot(1)*deltaT;
	
	for(i=2:(nMax-1))
		xDoubleDot(i) = (1/m)*( - Fapp(i) - ( k*(x(i)-(x(i-1)+l0))) - (b*xDot(i)) + k*(x(i+1)-(x(i)+l0)) );
		%Integration of Position and velocity;
		xDot(i) = xDot(i) + xDoubleDot(i)*deltaT;
		x(i) = x(i) + xDot(i)*deltaT;
		disp("inside loop");
	endfor

	%	xDoubleDot2 = (1/m)*(- Fapp2 - (k * (x2-(l0+x1))) - (b*xDot2));
	xDoubleDot(nMax) = (1/m) * ( -Fapp(nMax) - (k*(x(nMax) -(x(nMax-1) +l0)) - (b*xDot(nMax))  ));
	%Integration of Positionn and velocity;
	xDot(nMax) = xDot(nMax) + xDoubleDot(nMax)*deltaT;
	x(nMax) = x(nMax) + xDot(nMax)*deltaT;

%	disp("Booah");
	
	t= t+deltaT;

%	for ( i=1:nMax)
%		store(i,:)= [store(i,:); x(i)];
%	endfor

	%store = [: ; x(:)];
	store1 = [store1 ; x(1)];
	store2 = [store2 ; x(2)];
%	store2= [store2;x2];
%	plot (x(1),t,'r0'); hold on;
endwhile

%plot (store1, 'r');  hold on; 
%plot(store2 , 'b'); 
plot(store1,'r' );  hold on;
plot(store2,'b');